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Abstract. Quantum fluctuations of Bose-Einstein condensates trapped 
in disordered lattices are studied by inhomogeneous Bogoliubov theory. 
Weak-disorder perturbation theory is applied to compute the elastic 
scattering rate as well as the renormalized speed of sound in lattices of 
arbitrary dimensionality. Furthermore, analytical results for the con- 
densate depletion are presented, which are in good agreement with 
numerical data. 

1 Introduction 

In contrast to dense quantum fluids like liquid Helium, dilute gases of ultracold bosons 
can be condensed into rather pure Bose-Einstein condensate. But due to interactions, 
there is always a fraction of non-condensed particles, even at zero temperature [1-3]. 
This fraction can be enlarged by loading the atoms into the deep wells of optical 
lattices [4], until finally the celebrated quantum phase transition from superfluid 
to Mott insulator is reached [5—7]. A somewhat different quantum phase transition 
from the superfluid to the Bose-glass is driven by applying instead a potential that 
varies randomly in space [8—11]. Generally, the competition between kinetic energy, 
interaction and disorder makes this so-called dirty boson problem very rich in its 
phenomenology, and different parameter regimes require different approaches [12-14]. 

Here, we investigate some disorder-related properties of the deep superfluid phase 
preponderant at high filling in the Bose-Hubbard model, for which Bogoliubov theory 
provides the right framework [1, 15]. Bogoliubov theory has been used extensively 
for disordered Bose gases [16-18]. Here, we transfer our inhomogeneous Bogoliubov 
theory [19] from the continuum onto the lattice, which amounts mostly (but not 
entirely) to replacing the free-space dispersion by the lattice dispersion. The following 
sections 2 and 3 describe the model and its Bogoliubov theory in greater detail. As 
a first application, we compute the effective dispersion relation with scattering rates 
and speed-of-sound corrections in Sec. 4. 

It is an essential feature of our approach that quantum fluctuations are counted 
from the underlying condensate that is itself already deformed by the external poten- 
tial. In such a setting, it becomes quite a challenge to distinguish the condensed from 
the non-condensed component [20-23]. Recent Quantum Monte Carlo studies have 
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addressed this question both for weak lattices [24] and harmonic traps [ ]. In Sec. 5, 
we use the approach described in [26] to determine the disorder-induced quantum 
depletion in the lattice. We reach excellent agreement with numerical data obtained 
by Singh and Rokhsar [27] within an explicit Bogoliubov diagonalization, and thus 
validate our analytical formulation. Section 6 concludes. 



2 Bosons on a random lattice 

2.1 Disordered Bose-Hubbard model 

We consider bosonic atoms destroyed (created) by field operators c x (cj.) at the sites 
x = J2 < j=i n j a j> n j G Z, of a simple Bravais lattice with basis vectors a,, 1 < 
j < d, fulfilling the canonical commutation relations fecC^/l = S xx i, [c^cv] = 0, 
and [cjg, = 0. The system is governed by the grand-canonical single-band Bose- 
Hubbard Hamiltonian 

E = J2\ d x T xx>c x > + {V x - n)c x c x + ^c x c x c x c x \. (1) 

x ^ x ' ' 

The hopping matrix T xx > = - J2j=i Jj{Sx,x'+ aj + 5x,x'- aj - 25 x>x >) contains the 
tunneling amplitudes Jj along direction j between nearest-neighbor Wannier orbitals 
in the lowest band. To ease notations, we consider in the following a simple cubic 
lattice with edges of length Lj = Njd and a total number of sites N = YljJ\fj\ 
the generalization to other geometries is straightforward. The chemical potential [i 
controls the particle number, and U > is the onsite repulsion. Eq. (1) describes 
the simplest version of the disordered Bose-Hubbard model, where only the on-site 
energy V x is random [8,12,27]. More generally, randomness can also occur in hopping 
and interaction [21,22]. 

Without the external potential, the system is translation invariant, so it is appro- 
priate to use the Fourier representation Ck = 7V~ 3 ^2 X ex P (—ik ■ x) c x , with allowed 
wave vectors k in the first Brillouin zone (BZ) such that J2k ex p[*fe' ( x — x ')] — N&aias' ■ 
The kinetic energy operator becomes diagonal, and its matrix elements Tf-k' = &kk'£k 
define the clean lattice dispersion 

d 

4 = 2J^[l-cos(a%)]. (2) 

3=1 

In the continuous limit a — > 0, the hopping matrix turns into the continuous Lapla- 
cian, and 2Ja 2 —> ti 2 /m, such that e k — > h 2 k 2 /2m. In the remainder of this con- 
tribution, we will be concerned with the inverse mapping, namely to transfer the 
continuum results of [19,26] to the lattice case. 

2.2 Statistical properties of the disorder potential 

Hitherto, cold-atom experimentalists have explored three ways to realize random lat- 
tices: (i) Sampling a continuous potential V(r), such as optical speckle, at the lat- 
tice positions [21,22]; (ii) Trapping impurity atoms at or near the primary lattice 
sites [28,29]; (iii) Creating a quasi-random distribution by adding an incommensu- 
rate lattice to the primary lattice [30]. 
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Fig. 1. Schematic view of a disordered lattice hosting an extended, but deformed condensate. 
In the atomic limit J = 0, the condensate density (blue) would fill the wells up to (fi — V^/U. 
Finite hopping J (in this plot J = Un c ) couples adjacent sites and smoothens the profile. 
Sec. 3 and following discuss the quantum fluctuations 8h (and Sip), here cartooned in green. 



The random potential V x is characterized by its moments V Xl V X V X ', etc. Without 
loss of generality, we can set V x = by readjusting the zero of energy. Concerning the 
covariance, one is capable of realizing, in all of the above-mentioned cases, spatially 
uncorrelated disorder such that V X V X > — V 2 5 xx r, or 

VkV^ = V 2 ^- 1 S hk >, (3) 

for the Fourier components = Af^ 1 J2 X ex P ( — ^ ' x ) ^a>- ^ n the speckle case (i), 
this is valid if the spatial correlation length is not larger than the lattice constant; 
for the Bernoulli disorder (ii) [11,28], this is true by construction, and in a quasi- 
random system (iii) this requires the incommensurability to be large enough so that 
no repetitions occur in the finite size under study [31]. 

2.3 Condensate deformation on the mean-field level 

When the discreteness of the particle number at a given site is negligible, i.e., for 
high particle numbers per site, n ^> U/J, a mean-field description for the extended 
condensate is a good starting point (see the schematic representation in Fig. 1). 

Replacing the operators Cx by their ground-state expectation value <P X (which can 
be chosen real) and minimizing (1) yields the Gross-Pitaevskii equation, which we 
write in the Fourier representation: 

V-4>k = Efe^fc + ^2 V k-k'<l>k' + ^2 ^k-k'^k'-k'^k"- (4) 
k' k'k" 

If the external potential is weak, the homogeneous solution (f>^ = N^ 2 Sko is only 
slightly perturbed, and one can compute the deformed condensate amplitude 4>k — 
+ <p£' + . . . perturbatively in powers of V [12]. At fixed number of condensed 
particles, J2k l^k| 2 = ^c, also the chemical potential has to be expanded. As a re- 
sult, the condensate deformation reads <f>^ = —Nc Vk to first order in the reduced 
potential matrix element 



(5) 
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This expression matches Eq. (9) in [19] [or Eq. (15) in [12]], after substituting gn by 
Un c = UN C /Af, and reading er£ as the lattice dispersion (2). Also all higher-order 
terms agree with their respective expressions in [ 9], as do derived quantities like the 
inverse condensate amplitude and densities. 

In the limit J -C Un c of weak coupling between sites (which requires large occu- 
pation such that still Jn c 3> U), the bare kinetic energy e k in the denominator of the 
reduced potential (5) is always negligible compared to 2Un c . Thus, the condensate 
faithfully follows the potential, to all orders, resulting in a Thomas-Fermi profile and 
a screened disorder landscape [12]. In the strong-coupling limit J 3> Un c , kinetic and 
interaction energy become comparable at a characteristic wave vector equal to the 
inverse of the healing length, £ = a^f J J (Un c ) >• a. Whereas potential fluctuations on 
large scales 3> £ are still followed by the condensate, fluctuations on small scales 
k^ 1 -C £ cost too much kinetic energy, and result only in a smoothed imprint [32]. 

The deformed condensate consists of particles with the (ensemble-averaged) mo- 
mentum distribution [261 



nek = = N c (1 - TOw, + |V fc 2 | ; (6) 



the second equality holds to order V 2 in the disorder strength. Here, V% = J2k l^fel i s 
the fraction of condensed particles with k ^ 0. In the weak-coupling limit J <C Un c 
and for uncorrelated disorder (3), this fraction is V2 = v 2 /4, where v 2 = (V/Un c ) 2 
is the disorder variance in units of the chemical potential of the homogeneous lattice 
gas. 



3 Inhomogeneous Bogoliubov Hamiltonian on the lattice 

In this section, we set up the inhomogeneous Bogoliubov Hamiltonian for quantum 
fluctuations 6c x around the (deformed) mean-field condensate <P X on a disordered 
lattice. For the reasons given in [19], a well-defined theory ensues if one expresses the 
fluctuations via density and phase, 

Sti 

5c x =c x -<P x = + i<P x 5ip x . (7) 

Inserting (7) into (1), and collecting all second-order terms, 1 one arrives at a fluctu- 
ation Hamiltonian of the form 




Crossed terms have canceled, up to constants from commutators, and the coupling 
matrices between different Fourier modes are found to be 

S k k' = ^k-p^p-k* (e° - e° k _ p ), (9) 

c v 

2 

Rkk> = -jyT 2J 4>k-p<t>p-k> -q [(s p ~ n)6 q o + V q ] + 6Un c 6 kk > ■ (10) 
c p,q 

1 Higher-order terms describe interactions among the excitations that remain negligible 
for very low temperatures and dilute gases. 
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Fig. 2. (a) Equipotential lines of the clean and Bogoliubov lattice dispersions, (2) and (12), 
respectively, in the first BZ of the 2D lattice, (b) Plots of several quantities for J = 0.5t/n c 
along the closed path connecting F, K, and M in (a): Clean Bogoliubov dispersion eu, 
eq. (12), in units of 20Un c (full black line), together with the phonon dispersion Eh = hck 
for fc£ <C 1 (dashed black). Colored curves show the disorder corrections calculated in Sec. 4: 
relative scattering rate 7fc/efc (red); relative dispersion shift Ae^/Ek (blue), both in units of 
v 2 = V 2 /(Un c ) 2 . 



Here, 4>k = [n c /<P]k are Fourier components of the inverse condensate amplitude. 
Eqs. (9) and (10) are the lattice equivalents of (29) and (30) in [19]. There, we used 
the identity (j>V(j) = — </>V0 together with (4) to express Rkk 1 as a quadratic functional 
of 4>k without appearance of V q nor /i, which turns out impossible on the lattice. 

3.1 Homogeneous case 

In absence of the potential V x , the condensate is homogeneous, and the coupling 
matrices reduce to S^, = 2e1Skk' and R k °^, — 2[e" + 2Un c ]Skk>- The Hamiltonian 
becomes diagonal, = ^2 k EkT^k, after a canonical transformation to Bogoliubov 
quasiparticles ("bogolons") such that [7fc',7fc] = Skk'- This is achieved by 

(£)=M^%)< 

with parameter a,k = (efe/^fc) 1 / 2 in each fc-sector. The resulting Bogoliubov dispersion 

Sk = [el{el + 2Un c )] 1 ' 2 (12) 

is plotted in Fig. 2(b) along a representative path in the BZ of a 2D simple length 
£ = 2a/7r separates the sound-wave regime fc£ <§; 1 with linear dispersion from the 
particle-like excitation regime fc£ 3> 1 with parabolic dispersion. 



3.2 Bogoliubov scattering vertex 

For any given realization of the external potential V x , and thus a given realization of 
the coupling matrices (9) and (10), one can diagonalize the quadratic Hamiltonian (8). 
Analytical calculations for weak disorder, however, are best done in the bogolon basis 
(11) — with density (and phase) fluctuations counted from the deformed condensate. 
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The price for keeping a plane-wave basis is the appearance of an impurity-scattering 
term in the Hamiltonian, which can be treated using standard perturbation theory. 

We separate the homogeneous contributions from the coupling matrices by defin- 
ing RuOi = Rkk'—Rkk' > as we ^ as ^kO = ^fefe' — Suui ■ By the canonical transformation 
(11) to Nambu pseudo-spinors P k = (7L j- k )/y/2, the fluctuation Hamiltonian (8) 
is brought into impurity-scattering form [19]: 

h = J2 £ krlr k + flVkk>r k '- (13) 

k k,k' 

The bogolon scattering vertex 2 V kk > = {A k 1 )^ diag^S^J , R k V k i)A k } i s a function of 
the condensate ground-state configuration via (9) and (10), and can be expanded in 
powers of the bare potential using the perturbation theory of Sec. 2.3. Thus, our ap- 
proach goes beyond earlier Bogoliubov theories on disordered lattices [12,27] inasmuch 
our Hamiltonian is useful even for purely analytical calculations. To first order, 

R^l=2V k . k ,[e k _ k ,+e k + s° k ,], (14) 
S i k l^2V k . k> [e k _ k ,-e k -s k ,}. (15) 

If one is interested in results to second order at most, only disorder-averaged, diagonal 
terms are needed: 



S k l=2S kk ,^2\V q \ 2 [e°k- q -e k -s q ] (16) 
q 

R k l = 2<w E N 2 [4- q + $4 + 34] , (17) 

q 

Contrary to the continuous-limit case [19], where S kk ' of eq. (28) is proportional to 
^cfe-fc'j here its diagonal elements have no reason to vanish. 



4 Effective medium theory 

One of the first things one may want to know is how the disorder modifies the exci- 
tation spectrum. Clearly, randomness limits the life-time of plane-wave excitations, 
affects the speed of sound, changes the density of states, and causes localization 
[19, 3J 36]. All these effects can be assessed by applying diagrammatic perturbation 
theory to the Hamiltonian (13), for instance via a Green's function formalism. The 
Nambu-Green function of the clean Hamiltonian is Gok(z) — diag(Gofe(2), Go k (— z)), 
with Go k (z) = [z — £fc] _1 , and the full Green function obeys the recursive equation 
G = Go + GqVQ. The ensemble-averaged Green function G = Go + Go^G describes the 
propagation of bogolons through an effective medium, whose effect is to renormalize 
the dispersion: 

tw k = £fe + Eii(ftu, k) = s k +Ae k — 17*, (18) 

£11 designates the first block of the Nambu self-energy. To second order in the bare 
potential, one finds two types of terms, — V^GqV^ + V( 2 ), resulting from the 
perturbative expansion of the bogolon vertex. 

2 Eq. (31) in [19] contains a misprint with two factors of 1/2 too much; the following 
expressions are not affected. 
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Let us first discuss the elastic scattering rate 7&. To lowest order, known as the 
Born approximation, it is given by Fermi's Golden Rule via a momentum integral 
over the energy shell: 



^ fe = is 51 W a P R k ] p + a k la p ls S\ 2s ( £ k ~ Bp). (19) 
p 

Due to the rather complicated form of the energy shells, even for a simple cubic lattice, 
the integral has to be evaluated numerically. 3 In Fig. 2(b), the relative scattering rate 
Jk/£k in units of v 2 is shown across the BZ. Just as in the continuous limit, the 
scattering of low-energy excitations near the F point is suppressed, resulting in long- 
lived, soft modes. Strong scattering occurs near the symmetry point K at the BZ 
boundary. 

The real part of the self-energy in (18) shifts the energy of excitations. In Fig. 2(b), 
the shift Ae^/sk relative to the clean dispersion is plotted in units of v 2 . The shift As/~ 
is a principle- value integral over a kernel with a simple pole at the energy shell [19]. 
As shown in Fig. 2(a), the energy shell changes its topology at e = [4J(4J + Uric)} 1 / 2 , 
which causes jumps in the shift of the dispersion Ask- For the soft modes, Ask 
results in a changed speed of sound, c = c + Ac, which we observe to be negative in 
the 2D example of Fig. 2(b). In the limit J/Un c — £ 2 /a 2 3> 1, the lattice constant a 
drops out, and we recover the asymptotics for uncorrelated disorder in the continuum, 
Ac/c (x v 2 = v 2 (a/£,) d , where the proportionality factor is negative in ID, positive 
in 3D and vanishes in 2D [19]. The increase in 3D, Ac/c — 5u 2 /(48v / 2tt), coincides 
with the established result [17, 37, 38] for uncorrelated disorder. However, as noted 
in [19], an increase of c with disorder is untypical. Generically, one expects the disorder 
to reduce the speed of sound, and consequently to increase the low-energy density 
of states. And really, in the opposite, weak-coupling limit J <C Un c , the shift is 
negative and reduces to the value Ac/c = — 0.375u 2 , — 0.215v 2 , — 0.177i> 2 in dimensions 
1, 2, and 3, respectively. This is very similar to the hydrodynamic regime of the 
continuum [19,34], but with slightly different constants. 



5 Condensate depletion 

Due to the repulsive interaction between atoms, there is always a certain fraction 
of particles outside of the condensate. In a weakly interacting gas, this fraction can 
be estimated using Bogoliubov theory. Let us first discuss the homogeneous lattice 
case. Condensation occurs in the state with quasimomentum k = 0, and the non- 
condensed particles have a momentum distribution Sn^ = (Sc k Sck) = (flfe — a^ 1 ) 2 /4, 
as follows from (7) and (11) in the Bogoliubov vacuum at T = 0. The total density 
of non-condensed particles (number per lattice site) is 

5n^ =Af- i y ^ £k ~ £ p 2 . (20) 
k 4£k£ k 

In the continuous limit, J Un Cl it shows the well-known scaling proportional to 
{Un c /J) d l 2 = (a/£) d . In the weak coupling regime, J <C Un c , the asymptotics is 



3 We parametrize the integral by the polar (and azimuthal, in d — 3) angle, determine the 
modulus of p such that e v = Ek, and evaluate the integrand there. This has to be multiplied 
with the surface element and the Jacobian [9e p /9p„] _1 , where p„ is the component normal 
to the energy shell. For high energies, we take the corner M as origin. 
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V/Un c ,J/Un c 
(a) (b) 

Fig. 3. (a) 2D condensate fraction N c /N as function of the disorder strength v — V/Un c for 
different values of Un c / J. The data points are taken from Fig. 4 in [27], the result of an exact 
Bogoliubov diagonalization. The lines show our results for the condensate fraction resulting 
from the combined clean and potential depletion Sn = <5n' ' + 8nS 2 \ Eqs. (20) and (21). For 
not too strong disorder v = V/Un c < 1.5 the agreement is excellent, (b) Disorder-induced 
depletion, relative to the clean depletion, in dimensions d = 1, 2, 3. The dashed lines show the 
large- J asymptotics Sn^ /SrS ^ ~ fidV 2 = v 2 (Un c / J) d ^ 2 , corresponding to the uncorrelated 
limit of Ref. [26] . In one dimension, an infrared cutoff a = fcm£ is needed. 



SnW w g dy /Un c /J, where g 3 = 0.161 and g 2 = 0.227. In one dimension, an in- 
frared cutoff is needed, because there is no true Bose-Einstein condensation, and 
Sn^ grows logarithmically with the cutoff. The condensate fraction, or relative num- 
ber of particles in the condensate, then is N c (0)/N = 1 — 6nS°> jn. Using (20) in the 
thermodynamic limit, we find excellent agreement with the clean condensate fraction 
shown in Fig. 4 of Singh and Rokhsar [27], and infer n = 0.33 particles per site for 
their data. Even though the density is quite low, and the system parameters do not 
verify the prima facie condition U <C Jn of validity for mean-field theory, at this 
fractional occupation, the clean system has a superfluid ground state. Nonetheless, 
we cannot expect Bogoliubov theory to yield accurate predictions as the condensate 
fraction falls below 50% [4, 15]. 

Next, we study the excitation momentum distribution in the presence of disorder. 
We can express 5nk, still via (7) and (11), in terms of the condensate profile <£ x 

and bogolon ground-state expectation values such as ('Jk'Jk')', f° r details see [26]. To 
second order in the disorder potential, the result can be written 

p 

~ (2) 

Now the kernel M^ p of Eq. (49) in [26] has to be evaluated using the lattice dispersion 
(2) and the envelope functions that follow from (9) and (10). Just as in the contin- 
uous case of [26], the momentum distribution of excitations induced by the disorder 
broadens the condensate momentum distribution, Eq. (6). 

The resulting, additional quantum depletion of the condensate density caused by 
the disorder potential, or potential depletion for short, is the sum over the correspond- 
ing momentum distribution, 8n^ = Af^ 1 ^2 k Sn^ . In Fig. 3(a), we show results for 
two dimensions and compare them to data from the exact Bogoliubov diagonaliza- 
tion of Ref. [ ] . Within our perturbation theory, we can only estimate the quadratic 
corrections of order V 2 , but this appears to be a rather good approximation even 







for strong disorder up to V ~ Un c . The good agreement validates our analytical 
approach to inhomogeneous Bogoliubov theory. 

In contrast to [27], our theory is analytical up to the evaluation of the integrals, and 
we do not need to average over realizations of disorder. This allows us to consider 
higher dimensions with little additional effort. In Fig. 3(b), the relative potential 
depletion is shown as function of J/Un c in dimensions d = 1, 2, 3. In all cases, the 
depletion takes its maximum around J rs 0.5Un c . The large- J behavior coincides with 
the limit a — > of the continuous case. In accordance with [261, we find SnW 1 5 
favj, with fa « 0.236 (slightly cutoff-dependent), fa ~ 0.135, and fa « 0.160. 

In the weak coupling regime, J <C Un c , 8n^ shows the same scaling proportional 
to ^Un c /J as 5rv-°\ such that the ratio 5n^ / v 2 Sn^ tends to a constant (« 0.016 
in 3D, and w 0.06 in 2D). In one dimension, however, the weak-coupling behavior is 
cutoff-dependent: as J/Un c gets smaller than the cutoff a, the ratio diminishes and it 
is not possible to obtain a well-defined limit, as becomes evident from the two curves 
with different cutoffs shown in Fig. 3(b). 

To summarize this section, we have shown that the disorder potential does indeed 
increase the condensate depletion, essentially because it deforms the condensate and 
creates high-density regions, where the interaction-induced depletion is enhanced. 
The total depletion can be written [26] 

Sn = Sn^ [l + v 2 A(J/Un c )] , (22) 

where A(J/Un c ), plotted in 3(b), is always smaller than unity. Thus, the potential 
depletion for not-too-strong disorder remains a relatively small correction, comparable 
to the homogeneous depletion, up to V ~ Un c at least. For a large mismatch between 
J and Un c , the condensate becomes even more resilient to the disorder by virtue of 
smoothing (for J Un c ) and screening (for J -C Un c ). 



6 Conclusions 

We have studied the superfluid phase of the disordered Bose-Hubbard model by Bo- 
goliubov theory, transferring our previous continuum formulation to the lattice case. 
The key point is to describe the quantum fluctuations around a condensate that is 
deformed by the random potential. This approach permits purely analytical calcu- 
lations, at least to second order in the disorder strength. We have determined the 
renormalized dispersion relation, and find strong scattering at the K point (3D: X 
and M points), where the free dispersion relation has a saddle point. Otherwise, the 
excitations stay gapless, and soft modes are long lived, with well-defined corrections to 
the speed of sound. Furthermore, we have calculated by how much the non-condensed 
fraction increases due to the disorder. This additional depletion is found to agree with 
earlier, numerical results by Singh and Rokhsar. We determine its dependence on the 
system parameters and find that the disorder-induced depletion relative to the clean 
depletion is greatest when the hopping matches the chemical potential, which should 
provide useful information for further experimental studies. 
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